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Abstract 

We discuss the application of the static path plus random phase approximation (SPA+RPA) 
and the ensuing mean field+RPA treatment to the evaluation of entanglement in composite quan- 
tum systems at finite temperature. These methods involve just local diagonalizations and the 
determination of the generalized collective vibrational frequencies. As illustration, we evaluate the 
pairwise entanglement in a fully connected XXZ chain of n spins at finite temperature in a trans- 
verse magnetic field h. It is shown that already the mean field+RPA provides an accurate analytic 
description of the concurrence below the mean field critical region (|6| < 6c), exact for large n, 
whereas the full SPA+RPA is able to improve results for finite systems in the critical region. It is 
proved as well that for T > weak entanglement also arises when the ground state is separable 
> be), with the limit temperature for pairwise entanglement exhibiting quite distinct regimes 
for |6| <bc and |6| > be- 
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I. INTRODUCTION 



It is now well recoErnized that quantum entanerlement plays an essential role in both 
,.ant.. .fo..a.ion science g wh.e . i. co.ide.ed a a. well . . .a„.- 

body and condensed matter physics, where it provides a new perspective for understanding 
quantum correlations and critical phenomena |2|45[. Entanglement denotes those correlations 
with no classical analogue that can be exhibited by composite quantum systems, which 
constitute, for instance, the key ingredient in quantum teleportation [6]. A pure state of a 
composite system is entangled if it is not a product state, while a mixed state of such system 
is entangled when it cannot be written as a convex combination of product states 3]- 

Thermal entanglement [2], Isij]^ denotes that of mixed states of the form p(T) oc 
exp[— where H is the system Hamiltonian and /3 = 1/kT the inverse temperature. 
A complete characterization of thermal entanglement in many component systems is diffi- 
cult, since, to begin with, there is no simple necessary and sufficient computable criterion 
for determining if a general mixed state is entangled 13|. Besides, these systems exhibit 
entanglement at different levels, i.e., between any pair or set of subsystems, starting from 
that between elementary constituents i,j and ending in that of global partitions [14] (which 
for T > can no longer be measured through the entropy of a subsystem). Finally, a basic 
difficulty is the accurate evaluation of p(T) and the ensuing reduced densities pij. Standard 
methods like the mean field approximation (MFA), which may provide a correct basic de- 
scription of thermodynamic observables in some systems, are not suitable for the evaluation 
of entanglement since they are based on separable (non-entangled) trial densities. In small 
finite systems fluctuations of the order parameters become important [15] and the MFA 
is to be replaced at least with some average over different mean fleld densities, but such 
an approach will still fail to describe entanglement as it is essentially based on a convex 
combination of product densities. 

The principal goal of this work is to show the applicability of the static path plus random 
phase approximation (SPA-I-RPA) 16l-ll9| to the determination of thermal entanglement. 
The approach is derived from the path integral representation of the partition function 
obtained with the Hubbard- Stratonovich transformation 



20|, and has been applied to the 



description of 
physics [I6KI9 



Dasic observables in diverse fermionic models of nuclear and condensed matter 



21 



, l22| . It takes into account both the large amplitude static fluctuations 
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(SPA) of the mean field order parameters, essential in critical regions of finite systems, 
together with small amplitude quantum fluctuations (RPA), which may account for most 



quantum effects if T is not too low and wi 
fully consistent MFA+RPA approach 



1 be responsible for entanglement. It also provides a 
19| , obtained through the saddle point approximation 
to the full treatment. Here we will formulate the method for a system of n distinguishable 
constituents, where it involves in principle just local diagonalizations. 

We will employ the formalism to evaluate the thermal pairwise entanglement in a fully 
connected XX Z chain of n qubits or spins in the presence of a uniform transverse magnetic 
field h. Spin chains constitute an attractive scalable qubit representation for exploring 



and implementing quantum information processes [23|-|25| and can be realized in diverse 
physical systemSj_mcluding those based on quantum dots electron-spins {2^ and Josephson 
junction arrays [27], where the effective model includes coupling between any two spins. 
Fully connected symmetric spin models (simplex) have also intrinsic interest, providing a 
solvable scenario for examining entanglement in systems undergoing phase transitions. In 
particular, entanglement properties of the fully connected XX and XY model at T = were 



thoroughly analyzed in 28H30|. We will show that the XXZ model exhibits an interesting 



non trivial behavior at finite temperature, whose main features can be correctly described 
by the SPA+RPA for moderate finite n and even by the MFA+RPA below the critical 
region, the latter providing an analytic description which becomes exact for large n. The 
formalism is described in section II while application to the model is discussed in III. Finally, 
conclusions are drawn in IV. 



II. FORMALISM 

We will consider a composite system described by a Hamiltonian of the form 

H = - lY.^Q")" ^ (1) 

u 

where Q'^ are linear combinations of local operators, i.e., = Yli^i^ Q'^ ~ TliiQ^i^ 
with QJ' acting just on subsystem i (QJ' = Ji (g) . . . (g) QJ' (g) . . . O /„, with [Q^, Q'^] = 
if i 7^ j). In a spin chain Q'^ could stand, for instance, for total spin operators or general 
linear combinations Yli'^i^i '^^ individual spins s^. Any quadratic interaction between 
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subsystems, 

y=-lll O^vZ'Of , (2) 

where denote local operators, can be written in the diagonal form ([1]) (non-unique), after 
completing squares or diagonalizing the matrix Vi^ju' = v^,, with Q'^ suitable hnear combi- 
nations of the Oj^. We may assume > in without loss of generality if antihermitian 
operators Q'^ are allowed {vp{Q^Y ^ —Vv{iQ'^Y). In what follows we will consider finite 
Hilbert space dimension. 

The Hubbard-Stratonovich transformation allows then to express the partition function 
Z = Trexp[— as the path integral j2o| 

Z = I D[x] Trf exp{-^V[^ ^ + h[xir)]} , (3) 

h{x) = J2Hx)^ K{x)=H^-J2^-^Q"^ (4) 

i V 

where T denotes time ordering and the normalization JZ)[x] exp[— Jj|^(ir^j^x^(r)/2f,^] = 1 
is assumed. The integrand in (j3]) is essentially the trace of the imaginary time evolution 
operator U\x\ associated with the path x(r) and the linearized Hamiltonian /i[x(r)], and is 
here a product operator Ui\x\, not necessarily positive. Eq. ([3]) can be evaluated by means 
of a Fourier expansion 

x,{r) = x,\^xle'^-\ 00^ = 27171/13, (5) 

where the static coefficients, representing the time average (xi.(r))[o,/3], with 

D[x] oc UudxuUn^odxt 



In the SPA+RPA [l6|-ll9| (to be denoted for brevity as CSPA (correlated SPA)), the 
integrals over the static coefficients Xj, are fully preserved, while those over x", n 7^ 0, 
are evaluated in the saddle point approximation for each value of the x,^. The aim is to 
take into account large amplitude static fluctuations, which are particularly relevant in the 
transitional regions of finite systems, together with small amplitude quantum fluctuations, 
which should in principle account for most quantum effects if the temperature is not too 



low. The flnal result can be expressed as [19 1 

/oo 
e-^^^-i/'^''Zix)CnPAix)dix), (6) 
-00 



^CSPA 
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where d{x) = Yl^ ^/pJJ2^^uJ)dxy and 

Z{x) = Tr exp[— = JJ^ tr exp[— , (7) 

i 

oo 

CRPA(a;) = ]^Det[5i.,,/ + (8) 

n=l 

Ruu'ix,uj) = > ■ -. (9) 

with tr the local trace, \ki) the running local eigenstates {hi{x)\ki) = e^Jfej)) and p^. = 
g-/3efe, ^^j,g-/3/i.i(a:)_ gqg^ ([Qj)_([9j) involve just /oca/ diagonalizations. Eq. (|8]) is the RPA correc- 
tion, fundamental in the present context, which can be further expressed as |l7. Ifsj] 

where a = (/ej, k'^) runs over all pairs ki ^ k[ (a > indicating ki > /c-), Aq = e^- — and 
are the running RPA energies, determined as the roots of the equation 

Y)ei[5^yi + VyRuyi {x, oo)] = . (11) 

They come in pairs of opposite sign and can also be obtained as the eigenvalues of the matrix 

Aaa' (x) = XJaa' + Pa"^ V^Q'i^Q'"^, , (12) 

V 

where = pfc- — Vk\-, Q^a = (kilQiWi)- Eq. can be applied provided CupAix) > 0, which 
implies u"^ + uf > \/a,x. Since the lowest RPA energies Ua may become imaginary or 
complex for x away from the stable mean field solution (see below), the previous condition 
sets up a breakdown temperature T*, normally low, such that Eq. is applicable for 
T > T*. Setting Cb.pa{x) = 1 in ([6]) leads to the plain SPA which, although significantly 
improving the MFA in critical regions, is unable to describe entanglement, as it averages 
correspond essentially to those of a convex combination of separable densities (for h{x) 
hermitian) . 

MFA+RPA. Away from critical regions, we may also apply the saddle point approxima- 
tion to the static variables x^. This leads to the MFA-I-RPA (to be denoted as CMFA), 
given by pj] 



^CMFA 



e-(^^^^^'^/^^^Z{x)Co{x)C^PA{x) , (13) 



where x is the value which minimizes the "separable" free energy J-'(x) = ^yXl/2vij — 
TlnZ(x) and is determined by the self-consistent "Hartree" equations 

= v.(gOx, (14) 

with {Q^)x = J2i kPki{f^i\Qi\ki)- Co{x) accounts for the small amplitude static fluctuations 
and is given by 

= Bet[5,,, + v,{K,ix,0) -J2{k^mk^)^T'^\ (15) 

i,k ^ 

with dpkjdxy, = (3pk, T.k'{K\Qi'\K)i^kk' - Pk')- 

Away from critical points Eq. (1131) can be employed right up to T — 0. Note, however, 
that if the solution of exhibits a continuous degeneracy (due to a continuous symmetry 
violation by h{x)) the previous approach should be applied just to the intrinsic variables (see 



sec. III). In this case the lowest RPA energy vanishes at the mean field solution [18|, ll9| but 
Eq. f|T3l) is still applicable, as C^^p^{x), Eq. ([8]), remains finite for Ua 0. Omitting Crpa(3^) 
and Cq{x) in (fT3|) leads to the plain MFA, which corresponds to a separable (product) density. 

We may then employ Eqs. (E]) or ( |T3i) to calculate the two-site averages {0^0^') = 
{2/ f3)d\nZ/dvl^^, required to evaluate the reduced density pij and hence a certain monotone 
or measure of the entanglement between subsystems i and j. If not present in the original 
interaction, we may in principle add the necessary terms in V and set at the end v^^^^, = 0. 

III. APPLICATION 

A. Fully connected XXZ Model 

We will consider n qubits or spins coupled through a full range XXZ type interaction in 
the presence of a transverse magnetic field b. The Hamiltonian reads 

n n 

H = bY,st-VY,[s^s^^+sysy + {l-^)s^,s^^] (16a) 

i=l i^j 

= bS, - V[Sl + si + (1 - ^)Sl] + Eo , (16b) 

where denotes the spin at site i (considered dimensionless) , S = ^^e total spin, 

7 the anisotropy and Eq = nV{3 — 7)/4. It is apparent that H commutes with and 



S"^ — si + Sy + S"^, its eigenvalues being 

EsM = bM- V[S{S + 1) - 7M2] + Eo , (17) 

where M — —S, . . . ,S and S — S, . . . , n/2, with S — {^) for n even (odd). The ensuing 
partition function is 

n/2 5 

Z = TTexp[-/3H] = J]F(5) J] e-^^^^ , (18) 

where Y{S) = („/2-s) ~ (n/2-5-1)' with 1^(|) = 1, is the muhiphcity of states with total 
spin S and 5*^ — M, such that ^^5=5 i^('S')(25' + 1) — 2". In what follows we will write 

V^v/n, (19) 

such that all intensive energies Esn/n remain finite for n — )> 00 and finite v. 

We will analyze here the attractive case v > Q (and 7 < 1), where the ground state has 
maximum spin 5" = n/2 V 6, 7. If 7 < 0, the ground state will be fully aligned (|M| = n/2) 
V 6 7^ and no ground state entanglement will arise, whereas if 7 > 0, the ground state will 
exhibit as h increases n transitions M — >■ M — 1 at 

hM^lv{l-2M)/n, (20) 

where Esm — Es,m-i, becoming fully aligned for 

\b\>bc = jv{l-l/n). (21) 

Thus, be is the limit field for entanglement atT — 0, as all ground states with S — n/2 and 
\M\ < n/2 are entangled (see below). 



B. Exact concurrence 



We will examine here the entanglement of a pair of spins (i,j), which is determined 
by the reduced two-qubit density p2 = Pij — Tr n-{ij}p{T). In the present system p{T) is 
completely symmetric and p2 will obviously be identical for all pairs i ^ j- In the standard 
basis of sf , s| eigenstates, it will have the form 

/ p+ \ 

p a 

P2= , (22) 

a p 

^ p_ / 
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where p+ + 2p + j9_ = 1 and 

/^i \l2 =^ '^«/'^2 =^ '^J^'/ n(n-l) ^ 4 =^ n ' 

" V'^^ ■^J / n{n-l) ' 

with (O) = TrpO the thermal average and = ± isf. Hence, p2 is here completely 
determined by the three collective averages (5*2), {S^) and (5*^), which can be directly derived 
from Eq. ( ITSl) as (we set k = 1) 

We may equivalently use (5*^) = {nT / v)d In Z/d'j. 



As a measure of pairwise entanglement we will employ the concurrence C 311], which for 



a general two component system can be defined as the minimum, over all representations 



Pi3 = E.9-l^-)(*-l, of E.^-C'd^-)), with C(|^^)) = v/2[l-tr(pn'] the square root of 
the linear entropy of any of the subsystems [32]. The entanglement of formation |33| is 
similarly defined but with C(|\E'^)) replaced by the standard entropy — TrpJ' log2 pj'. 

For a two qubit system, C can be explicitly computed as [311] C = [2 A — tri?, 0]+, where 
[u]^ = |(u + |m|)/2 and A is the largest eigenvalue of = [■^/p2P2^/P2Y^'^ , with p2 = 
4"^ s^j P2S^j the spin flipped density. The entanglement of formation becomes then just an 



increasing function of C (given by = — X]i/=± li' I0S2 with g± = (1 ± y/l — C^)/2), with 
ii^ = C = (1) for a separable (maximally entangled) pair. 
In the present system we then obtain 

C = 2[|a|-v^]+ (24) 



Y - (5,)2]+ (25) 



n n — 1 V n — 1 

so that P2 will be entangled if and only if |q;| > ^/p+P- , a condition which directly follows 



from Peres criterio n [34jl . In (!25|) 2/ri is the maximum value that can be attained by C in 
symmetric systems jsS], reached here for S = n/2 and M = ±(r;,/2 — 1) (in which case \SM) 
is an ly-state). 

T = behavior. Let us first briefly discuss the concurrence in the T — )■ limit, where 
and Sz approach sharp values S{S + 1) and M, with 5* = n/2. Eq. ( I25l) becomes then 
almost constant except for \M\ close to n/2 — 1, leading, up to 0((n — 1)^^), to 

C ~ ^ + ^ (26) 

n - 1 1 - 4m2 (n - 1)2 ' ^ ^ 
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FIG. 1: (Color online) Concurrence C (multiplied by n) at low temperatures as a function of the 
magnetic field b for n = 20 spins coupled through the Hamiltonian (jl6p for 7 = 1. The curves 
depict exact and CMFA results for T/v = 0.005 (top panel), where C reaches its maximum value 
2/n for b ~ b^ and T/v = l/2n = 0.025 (bottom panel), where the peak at 6 ~ 6c is no longer 
prominent (Eq. (j4ip ). The T — )• behavior for any 7 > is identical except for the rescaling 
V — > 7f . 

for m = M/n ^1/2. C increases stepwise from l/(?2— 1) for M = to 2/n for \M\ = n/2 — 1, 
vanishing for \M\ = n/2 28|. The ensuing behavior of C for n = 20 is depicted in Fig. [TJ 
The dips occur at the field values fl20|) where the levels cross, in which case Eq. f l25|) leads to 
a strictly constant lower value C = 1/n due to the fluctuation (S*^) — {Sz^ = 1/4 at these 
points. 

Thermal behavior. The concurrence (1251) vanishes in symmetric states with fixed S and M 
if S* < n/2, with the only exception of the case \M\ = S = n/2 — 1, where C = 2/{n{n — 1)) 
(and a < 0). Hence, for \b\ < be we may expect a monotonous decrease of C with increasing 
temperature, as the essential contribution will come from the states with S = n/2. The 
behavior for \b\ < be will be discussed in detail in the next subsection. 

Nevertheless, for T > a weak pairwise entanglement also arises for \b\ > be, i.e., when 
the ground state is fully separable, up to a limit temperature Tl that becomes constant for 
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arge b. The behavior is thus similar to that arising with nearest neighbor XX couphng 
3 (and in agreement with the persistence of global entanglement for large fields in XX Z 

models [ij]), although here Ti will decrease as n ^ for large n with the scaling fll9p . To 
prove this result, we set 6 > and note that for b — be ^ T , we may just keep in Z states 
with zero, one and two spins up (M = — n/2 + 0, 1, 2) for evaluating C in the lowest non-zero 
order {0{e~'^'')). This leads to 



C ^ '^^I^l^ii _ ^-Pv _ l^^-Mn^^ ^ (27) 

n \ n — 1 

r, = l-{n- l)e-^" + ln{n - 3)e-2/3-(i-i/") . (28) 

The field dependence in this limit is thus reduced to an exponential decay, with the limit 
temperature b-independent and determined by the root of the bracket in fl^7|) (always positive 
for T — )■ if 7 > 0). For large n, C is positive just for low T oc and we may accurately 
neglect e~^^ and set 77 ^ 1 in fl27|) (in which case it is just the result from the S = n/2 
multiplet). This yields 

2'jv 2'jv 
'^'-'^ n\n[2n/{n-l)] "^^^2 + 1' ^ ' ^^^^ 

The maximum value reached by C in this region (attained close to Tl) is very small (oc 

^-2^-n{ln2)ib-yv)/2y 



C. CSPA and CMFA results for the XX case 

We start by describing the XX case (7 = 1 in ( fT6l) ). In the representation ( 116bl) . the 
CSPA, Eq. (|6]), will lead to a two-dimensional integral over variables {x,y) = r(cos0, sin0) 
associated with the linearized Hamiltonian h{x,y) = bS^ — xS^ — ySy + Eq/u. Since 
[H,Sz] = 0, both Z{x,y) and C^pA{x,y) will be independent of the orientation 0, and 
the final expression can be written as 

^CSPA = ^ /"Trfre-"^^'/^^Z(A)CRPA(A, u) , (30) 



where A = y/lP + r^ is the energy gap determined by h{x, y) and 



Z(A) =e-^^°(2cosh^)'^, (31) 
a;sinh(/3A/2) 

Asinh(/3c./2)' (^2) 



00= \/{X-vtanhf){X-v^tanhf). (33) 
10 



There is here a single collective RPA energy u. For temperatures lower than the mean 
field critical temperature Tc (see below), Eq. fl33|) becomes imaginary for r in an interval 
just below the stationary point (where w = 0), leading to the CSPA breakdown when 
cj^ < — 47r^T^. This is first satisfied at 6 = and r ^ f/2, where u ^ leading to a 

breakdown temperature T* v / 4:7^ {h = Q) . T* decreases as b increases, vanishing for b > v. 
CMFA. The mean field equations (HM reduce here to 

r = w— tanh^, (34) 
A 2 

and determine the minimum of the "Hartree" potential J'(r) = nrV(4f) -Tln2'(A). We 
then need to distinguish between two regimes: 
a) For |6| < v and T < T^., where 

^.= H/ln[^^, i\b\<v) (35) 

the minimum of J^{r) occurs at r > 0. This solution of f p4|) breaks the rotational symmetry 
around the z axes and is hence continuously degenerate (0 undetermined). In this case 
Eq. flM|) implies that A is the root of A = i; tanh(/3A/2), being hence b- independent, the 
constraint A > |6| leading to the critical temperature (!35|) (which is a decreasing function of 
\b\, vanishing for |6| — )■ v and approaching v/2 for b 0). The gaussian approximation (fT3|l 
in the "intrinsic" variable r leads then to 



^CMFA = e-i^(^^-^^)Z(A)sinh f , (36) 

X = |/3t;/cosh2f = iMl-S), (37) 

where the first two factors in (1361) represent the MFA result and the rest the RPA and SPA 
corrections. Note that u vanishes at this solution, in agreement with the broken continuous 
symmetry, but the RPA correction (|32|) remains finite (and essential) for a; — t- 0, with 
Crpa(A,u;) ^ sinh(/3A/2)/(/3A/2). 

It is apparent from fl30|) and fl3Tl) that in this region the approximation fl36|) will become 
increasingly accurate as n increases (the r fiuctuation decreasing as n~^), approaching the 
exact result for n — > oo. 

b) For |6| > f or T > T^, the minimum occurs at r = (normal solution of fl34p ). Direct 

application of Eq. (fT3|l in the original variables (x, y) leads then to 

_ sinh(/3&/2) 
^CMFA - ^(^)3i^h(/3a;/2)' ^^^^ 

11 



where u = b — t'tanh(/36/2). 

Evaluation of the concurrence. Let us first examine the CMFA concurrence for |6| < v. 
Both {Sz) and (S*^) (Eq. (12^ ) will be determined just by the MFA contribution in as 
the rest is b independent, and given by 

{S.) = -n^, {Sl) = {Sz)' + n^, {T<T^). (39) 

Hence, in this regime (Sz) is independent ofT (it is the value minimizing b{Sz) + v{Sz)'^ /n) 
while the fluctuation (S^) — (•S'^)^ increases linearly with T, reflecting a gaussian distribution 
p{M) oc e~'^^(*^~^'^^^)^/". In contrast, (S"^) is affected by all terms in ( 136|1 and given by 

being 6 independent. The first term is the Hartree part. For T— J-O, A/f— i-l while % — > 0, 
so that fHOj) approaches the right limit ^(^ + 1) owing to the RPA correction. 

The CMFA concurrence is then obtained replacing Eqs. ( !39l) -(l40l) in ( l25l) . As seen in 
Figs. (II])-(|2]), CMFA results turn out to be extremely accurate below the critical region, 
being undistinguishable from the exact ones if T is not too small. For T — )■ 0, CMFA 
actually leads to the exact result but with M replaced by the continuous variable ^nb/v, 
representing then the exact n — )■ oo limit. Accordingly, it does not reproduce the stepwise 
behavior arising for T — )■ and finite n, but remains close to the exact curve, correctly 
predicting the peak at 6 ~ 6c (top panel in Fig. [1]). 

For low T ^ Tc, thermal effects in the CMFA will arise just from the Sz fluctuation in 
fl5^ . as we may still set x = in (HUl) . As seen in the lower panel of Fig. [H CMFA correctly 
predicts the low temperature 

f = v/i2n) , (41) 

where the peak at 6 ~ 6c disappears. In fact, at T = T the CMFA concurrence has a 
strictly constant value C = 1/n for b < b^ while for T > T it starts to decrease with 
increasing field. We also note that for T < T the CMFA result is applicable just for 



b < b* = be — v\/l — T/T/n, becoming complex for b > b* and being maximum just at 



b = b*, where C = (l + yl — T/T)/n. For T > T it can however be applied right up to the 
limit field where C vanishes in the CMFA. 

As seen in Fig. [2l as T increases beyond T, the CMFA provides practically exact results 
for C even for n = 20 if |6| < since the concurrence in this region vanishes below 

12 
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FIG. 2: (Color online) Concurrence for n = 20 spins and 7 = 1, as a function of the magnetic 
field at different temperatures (top) and as a function of temperature at different fields (bottom). 
Exact, CMFA and CSPA results are depicted, that of CMFA vanishing for h/v = 1.1 (lower panel), 
where entanglement arises for T > 0. 

Tc. However, the CMFA accuracy decreases significantly if |6| is close to v. Moreover, for 
|6| > w the CMFA (Eq. fl55]) ) is not able to reproduce the exponentially small entanglement 
arising in this region. For large fields Eq. ( |38l) leads to an expression similar to ( 1271) . i.e.. 



always negative since it lacks the last exponential factor present in f l27|) . 

On the other hand, the full CSPA significantly improves CMFA results in the wide tran- 
sitional region around b ^ v arising for small n, as seen in Fig. [2j We may also appreciate 
the improvement over CMFA at field b = 0.9v, where the CMFA result is inaccurate for 
all T whereas the CSPA result is practically exact above the breakdown temperature, and 
also at 6 = l.lf , where the CMFA result vanishes while CSPA does predict the reentry of 
entanglement for T > 0, albeit above a certain onset temperature. Nonetheless, the CSPA 
result cannot reproduce the exponentially small entanglement arising for very large fields 
either, since it vanishes above a certain limit field larger than v. 



n Ln— 1 V V "~ 



with 77' — )• 1 for T <^ f , but the bracket is now 
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FIG. 3: (Color online) Top: Concurrence as a function of the magnetic field at T/v = 0.1 and 
different values of the number n of spins, for 7 = 1. Exact, CMFA and CSPA results practically 
overlap for n > 100, the concurrence vanishing for n > 8810. Bottom: Concurrence as function 
of temperature at fixed field b = O.Sv and increasing values of n. Exact and CMFA results are 
undistinguishable for n > 100. 

Results for large n. As n increases, the width of the transitional region diminishes and 
the CMFA prediction for |6| < f becomes increasingly accurate, being practically exact if 
n > 100, as seen in Fig. [31 An expansion of the CMFA concurrence up to Oin"^) leads to 



n 



:i-x)^ 



1 - 67^2 



1+ ) 



(42) 



where the first term contains the RPA correction and provides the only positive contribution. 
For T < Tc, Eq. ([MD leads to \/v^l- 26"^^ ^ 1 - 2e-^\ in which case x ^ 2/3^6"^^ and 
(142|) reduces, up to order x, to 



n 1 — o^/i;^ 



(43) 



Eq. (HSj) provides a simple yet accurate description of C for n > 100 if |6| <v. It implies an 
initial quadratic decrease with increasing field, and an initial n-independent linear decrease 
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T 1 1 1 1 1 1 — r"| 1 1 r 




b/v 



FIG. 4: (Color online) Limit temperature for entanglement as a function of the magnetic field for 
different values of n and 7 = 1. Exact, CMFA and CSPA (for n = 20 and n = 100) results are 
depicted, undistinguishable for n > 100. Tl decreases logarithmically with n ior b < v and as 1/n 
for b > V. The sparse dashed line indicates the mean field critical temperature Tc. 

with increasing T (Fig. [3]) for 2ne~^^ <^ 1, followed by a pronounced ra-dependent decrease 
arising from the exponential term in ( H3l) (which represents the effect from the S = n/2 — 1 
multiplet). Note also that at fixed T, entanglement will disappear for n > |e^^(l — IT jv) 
[n > 8810 in top panel of Fig. |3]). 

Eq. ( 1l3|) leads to a simple analytic expression for the limit field for entanglement &l(T), 

bL{T)^v\ l ^^^"^ ^ , (44) 

which is accurate for large n and T > (Eq. fl29|) ). as seen in Fig. HJ The inverse of (jH]) 
is the limit temperature TL{b), which is always lower than Tc (Eq. (|35i) ) and exhibits two 
regimes: For Tl(6) < V (2ne-^/^^(^) < 1), 



Tdb) ^v-\b\, (45) 
which applies for large in a narrow field interval just before |6| = v, whereas for large n, 

indicating a logarithmic decrease with n, in contrast with the decrease arising for |6| > v 
(Eq. fl2^ ). For very large n, this yields Ti(6) ~ v/\ia{2n), independent of b. We also note 
in Fig. m that for n = 20 and 100, the CSPA improves the CMFA prediction of T^ib) in 
the critical region, up to the field region where Tl(6) becomes close to the asymptotic value 
( 129|) . although it also leads to a lower onset temperature (lower sparse dotted lines). 
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D. CSPA and CMFA in the XXZ case 



In the general case (7 < 1) the CSPA leads to 



^CSPA 



/•oo /'OO ^ o 

/ rdr dz e-^(' +— ^Z(A)Crpa(A, (47) 

Jo J -00 



4 y 7rt>3(l — 7) 

where Z{X) and Crpa(A,Co') are given by Eqs. (ISI])-(IS2]) with A = ^{b-zy + r"^ and 



oj = -^/(A-t;tanh^)(A-i;(l -7§)tanhf ). (48) 

The mean field equations (IT^ become now 

r /3A b — z f3X 

r = V— tanh — , 2; = (7 — Ijt' — - — tanh — . (49) 
A 2 A 2 

In the symmetry-breaking phase (r > 0), feasible for 7 > 0, the solution for A is then 
identical with that for 7 = 1, i.e., A = f tanh(/3A/2), independent of b and 7, in which case 
( H9|) leads to 6 — 2; = 6/7, independent of T and v. This implies the rescaling 6 —t- 6/7 at 
the CMFA level. This phase is then feasible for |6| < 7^ and T < T^, where is given by 
Eq. ( l35|) with b — )• 6/7. The ensuing gaussian approximation to both z and r in ( l47l) leads 
to the CMFA partition function 

^cmfa(7,^^,^) = ZcMFAil,b/'j,v,T)/y/j , (50) 



where ZcMFA{'^,b,v,T) is the result 

From Eq. dSOD we obtain {S,) = -^nb/{-fv), and {S^,) - {S,)^ = \nT/{nv) [y -fv in 
fl39] and hence in f HTj) ). whereas (S*^) remains unchanged (Eq. fHOl) ). The ensuing results 
for C exhibit the same previous features, CMFA being accurate for \b\ < 7f , and the CSPA 
improving the latter in the transitional region |6| ^ 7^ . This can be seen in Fig. (|5]), whose 
upper panel depicts the quenching of the exact and approximate limit temperatures Tiip) 
for increasing 7. 

For large n and T <^Tc, the CMFA leads now to 



which generalizes Eq. ( l43l) and provides an accurate description for n > 100 if |&| < ^v. 
For 7 < 1 it implies a more pronounced initial linear decrease with increasing T, as seen in 
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FIG. 5: (Color online) Top: Limit temperature for entanglement as a function of the magnetic 
field b for different values of the anisotropy 7 in (jl6p and n = 100. Exact, CSPA and CMFA 
approximate results are depicted. Bottom: Concurrence C as a function of temperature at zero 
field for n = 20, 100 and 1000 and two different anisotropics. 

the bottom panel of Fig. [5l which for low 7 may persist up to the vanishing of C even for 
moderate sizes [n ^ 100 in Fig. [5]). Eq. fl5T]) leads to a limit field 



2T/M 



(52) 



1 - 2ne-l^'' ' 

which describes the CMFA results of Fig. ([5]). It implies 

T^(6)^7^'[1-&V(7^)'], 

for 2ne~'"^'^^^''^ ^ 1, a condition which may now apply V |6| < 'jv for moderate n if 7 is 
sufficiently low, whereas for sufficiently large n, 

2/7 



nib) 



■[1 



\n{2n)' (ln2n)2[l -62/(7t;)2]^ ' 

decreasing logarithmically with n and becoming independent of 7 and b (for |6| < 7^) for 
very large n. 
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IV. CONCLUSIONS 



We have shown the feasibihty of the CSPA approach for the determination of the pairwise 
entanglement in composite systems at finite temperature. The method is tractable, requiring 
just local diagonalizations and the evaluation of the generalized RPA energies, and unveils 
the crucial role played by the RPA correlations in the description of entanglement. It also 
leads to a consistent mean-field+RPA treatment (CMFA), which, as seen in the example 
considered, remains applicable and accurate in the presence of vanishing RPA energies, 
arising when continuous symmetries are broken at the mean field level. 

In the XXZ model considered, the CMFA provides an accurate analytic evaluation of 
the concurrence below the critical region well below even in relatively small systems, 
providing exact results for large n if |6| < 6c- The full CSPA allows to extend the accuracy to 
the critical region (6 ~ 6c) in finite systems, above a low breakdown temperature, predicting 
a reentry of entanglement for T > for fields above but not too far from he-, not detected 
at the CMFA level. Neither CMFA nor CSPA predict, however, the exponentially small 
entanglement arising for large fields at low non-zero temperatures. 

The present results also reveal the rich thermal entanglement properties of the fully 
connected XXZ model. We have shown by means of the CMFA that the limit temperature 
for pairwise entanglement decreases as f/(ln2n) for very large n if |6| < 6cj whereas for 
|6| > hf. it decreases as ^vjn^ becoming independent of h for large fields. CMFA also shows 
that the T = peak in the concurrence just before the transition to the aligned state at 
b — be disappears at a low temperature T ^ ^v/{2n). The extension of the present approach 
to more complex systems, including spin chains with general anisotropic interactions, is 
currently under investigation. 

The authors acknowledge support from CIC (RR) and CONICET (JMM and NC) of 
Argentina, respectively. 
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